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c3 Abstract 

We consider the problem of stochastic prediction and control in a time-dependent stochastic 

environment, such as the ocean, where escape from an almost invariant region occurs due to 

random fluctuations. We determine high-probability control-actuation sets by computing regions 
of uncertainty, almost invariant sets, and Lagrangian Coherent Structures. The combination of 
geometric and probabilistic methods allows us to design regions of control that provide an increase 
in loitering time while minimizing the amount of control actuation. We show how the loitering 
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^ time in almost invariant sets scales exponentially with respect to the control actuation, causing an 
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exponential increase in loitering times with only small changes in actuation force. The result is 
^ that the control actuation makes almost invariant sets more invariant. 
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Prediction and control of the motion of an object in time-dependent and 
stochastic environments is an important and fundamental problem in nonlinecir 
dynamical systems. One of the main goals of control is the design of a the- 
ory that can take unstable states and render them stable. For example, small 
perturbations at the base of an inverted pendulum will stabilize the inverted 
state. Noise poses a greater problem for deterministically controlled states, 
in that stochastic effects destabilize the states as well as their neighborhoods. 
Therefore, control theory of stochastic dynamical systems may be addressed by 
examining the change in stability of certain sets. 

We present a variety of geometric and probabilistic set-based methods that 
enable one to compute controllable sets. For a particle moving under the influ- 
ence of deterministic and stocheistic forces with no control, these sets determine 
regions which are unstable in the sense that the particle will leave the set after a 
sufficiently long time. Controls added to peirticle dynamics to increase the time 
to escape (or loitering time) have a strong dependence on the probabilistic and 
geometric set characteristics. The determination of controls and associated sets 
allow for an increase in the amount of time the particle can loiter in a particular 
region while minimizing the amount of control actuation. 

Our theoretical analysis shows how an increase in the strength of the control 
force leads to a decrease in the probability that an object will escape from the 
control region. In fact, we have found that small changes in the control actuation 
force have an exponential effect on the loitering time of the object. Additionally, 
we show how the exponential increase in escape times from the controlled sets 
is related to the problem of noise-induced escape from a potential well. 

I. INTRODUCTION 

Ocean circulation impacts weather, climate, marine fish and mammal populations, and 
contaminant transport, making ocean dynamics of great industrial, military and scientific 

interest. A major problem in the field of ocean dynamics involves forecasting, or predicting, a 
variety of physical quantities, including temperature, salinity and density. By fusing recently 
measured data with detailed flow models, it is possible to achieve improved prediction. 
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Therefore, forecasts occur with increased accuracj^. As a result, ocean surveillance may be 
improved by incorporating the continuous monitoring of a region of interest. 

Researchers have used surface drifters and submerged floats to acquire data for many 
years. More recently, sensing platforms such as autonomous underwater gliderP'Sl have 
been developed. The gliders can operate in both littoral (coastal) and deep-ocean regimes, 
and may be used for data acquisition, surveillance and reconnaissance. One drawback in 
the use of gliders involves their limited amount of total control actuation due to energy 
constraints, such as short battery life. For applications such as regional surveillance, energy 
constraints may be alleviated by taking advantage of the dynamical flow field structures and 
their respective body forces. 

Autonomous underwater gliders are subjected to drift due to hydrodynamic forces. This 
drift can be extremely complicated since the velocity fields found in the ocean are aperiodic 
in time and feature a complex spatial structure. Instead of constantly reacting to the 
drift (and thereby expending energy), one can minimize the glider's energy expenditure by 
taking advantage of the underlying structure found in geophysical fiows. However, in order 
to harness the ocean forces to minimize energy expenditure during control actuation, one 
needs to analyze the structures from the correct dynamical viewpoint. 

The potential for dynamical systems tools to shed light on complex, even turbulent, fiow 
fields such as ocean dynamics, has been understood for decades. Work in this area has 
intensified in the past decade with new focus placed on the improved Lagrangian perspec- 
tive that dynamical systems approaches may provide, especially for complex, aperiodic fiows 
in which the traditional Eulerian perspective on the fiow is unhelpful or even misleading. 
During boundary layer separation, for example, sheet-like fiow structures coincide with fiuid 
particles being ejected from the wall region, a technologically important feature that an 
Eulerian framework fails to capture in unsteady cases, but which has a clear Lagrangian sig- 
nature that can be identified using finite-time Lyapunov exponent^. In liquid jet breakup, 
a similar ejection process occurs during primary atomization, but a limited Lagrangian per- 
spective is provided by the liquid-gas interface, revealing liquid sheets and ligaments that 
precede droplet formatiorP. The development of these critical fiuid sheets and ligaments 
can be traced to unsteady, finite-amplitude global fiow structured. Dynamical systems 
tools thus provide a new approach to the characterization and control of these important 
fiows. Geophysical fiows offer another example where an Eulerian framework is ineffective 
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in the diagnosis of large Lagrangian structures and the measurement of transport. For pre- 
diction and control of particle dynamics in large surveillance regions of interest, Lagrangian 
structures of geophysical flows need to be characterized in both deterministic and stochastic 
settings. 

The field of Geophysical Fluid Dynamics (GFD) involves the study of large-scale fluid 
flows that occur naturally. GFD flows are, by nature, aperiodic and stochastic. The data sets 
describing them are usually finite-time and of low resolution. Established tools of dynamical 
systems have proven to be less effective in these cases: while providing some insight, they 
cannot provide realistic or detailed flow field data relevant to the trajectories of tracer 
particles. For this, dynamical systems tools can be applied to fluid flows in an alternative 
way, by interpreting the Eulerian flow field u(x, t) as a dynamical system x' = u(x, t) 
describing the trajectories of tracer particles. The phase space and real space are identical 
here and due to the incompressibility condition V ■ u = 0, the resulting dynamical system 
is conservative. This "chaotic advection" approach originated with AreP and demonstrated 
that even simple laminar two-dimensional (2D) periodic and three-dimensional (3D) steady 
flow^ could lead to complex, chaotic particle trajectories. Figure [l] illustrates an example of 
how particles in a periodic flow may exhibit unexpected trajectories. In this figure, a single- 
layer quasi-geostrophic beta plane is being driven by a bimodal wind-stress with a small 
amplitude periodic perturbation. Details of this ocean model can be found in Appendix |X} 

In the last two decades, this approach has led to new tools including the study of transport 
by coherent structured, lobe dynamic^, distinguished trajectories, and global bifurca- 
tion j^^ ' ^^ 4 Even though the transport controlling structures in GFD flows are inherently 
complicated and unsteady, their understanding is necessary to the design of glider controls. 
To overcome this obstacle we combine a set of dynamical systems tools that have proven 
effective in higher dimension^ and stochastic problemS^^. 

In this paper, we will consider a well-known driven double-gyre flow as an example to 
illustrate our prediction/control frameworld. This model can be thought of as a simplified 
version of the double-gyre shown in Fig. [T] which is a solution to a realistic quasi-geostrophic 
ocean model. It should be noted that our methods are general and may be applied to any 
flow of interest. The goal of our approach is the production of a complete picture of particle 
trajectories and tracer lingering times that enables one to design a control strategy that 
limits tracers from switching between gyres. 
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FIG. 1: (Color online) Streamfunction (color) and two particle trajectories within a single-layer 
quasi-geostrophic ocean model (see Appendix [A]) subject to low-amplitude periodic forcing whose 
mean gives a steady double-gyre/western boundary current flow solution. Initial points on the 
trajectory are identified with open circles. 

The techniques we will use to analyze the dynamical systems are based on both deter- 
ministic and stochastic analysis methods, and reveal different structures depending on the 
system under examination. In the deterministic case, the Lagrangian Coherent Structures 
reveal much about transport and information related to the basin boundaries. They also 
pinpoint regions of local sensitivity in the phase space of interest. Since basin boundaries are 
most sensitive to initial conditions, uncertainties in the data near the boundaries generate 
obstructions to predictability. Therefore, highly uncertain regions in phase space may be 
revealed by computing local probability densities of uncertain regions in the deterministic 
case, but using noisy initial data. Finally, in the time- dependent stochastic case, we can 
describe which sets contain very long term, albeit finite, trajectories. The sets are almost in- 
variant due to the stochastic forcing on the system, which causes random switching between 
the almost invariant sets. The tools we use here are based on the stochastic Frobenius- 
Perron operator theory. Once the full structure of almost invariant sets is identified along 
with regions of high uncertainty, control strategies may be designed to maintain long time 
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trajectories within a given region with minimal actuation. In Fig. [T| one can see that one 
particle's trajectory remains in the lower gyre for a long period of time, while the other par- 
ticle escapes from the lower gyre to the upper gyre. The tools we will develop and outline in 
this article will enable one to know if and when a control force must be actuated to prevent 
the particle from escaping. 

The layout of the paper is as follows. In Sec. |TT] we present the stochastic double-gyre 



system and examine the deterministic dynamical features. In Sec. Ill, we show how to use 



the finite-time Lyapunov exponents to describe transport, and we show in Sec. IV how to 



quantify uncertainty regions. We then turn to the stochastic system, and describe in Sec. [V] 



how to compute almost invariant sets. Section VI contains a discussion of our corral control 



strategy, and Sec. VII contains the conclusions and discussion. 



II. THE MODEL 

A simple model of the wind-driven double-gyre flow is provided by 

X = —nA sin(7r/(x, t)) cos(7r?/) — ax + r]i(t), (la) 

df 

y = 7TAcos{'Kf{x,t))sm{7Ty)— -ay + r]2{t), (lb) 

/(x, t) = e sm{ujt + ijj)x'^ + (1 — 2e sm{ut + ^/j))x. (Ic) 

When e = 0, the double- gyre flow is time-independent, while for e 7^ 0, the gyres undergo 



a periodic expansion and contraction in the x-direction. In Eqs. (la)-(lc), A approximately 
determines the amplitude of the velocity vectors, 1^/27? gives the oscillation frequency, e 
determines the amplitude of the left-right motion of the separatrix between the two gyres, 
is the phase, a determines the dissipation, and rii(t) describes a stochastic white noise 
with mean zero and standard deviation a = \/2D, for noise intensity D. This noise is 
characterized by the following first and second order statistics: {rii{t)) = 0, and {f]i{t)^j{t')) = 
2D6ij6{t — t') for i = 1,2. In the rest of the article we shall use the following parameter 
values: a = 0.005, A = 0.1, e = 0.15, and u = 27r/20. We consider the dynamics restricted 
to the domain {{x,y) | < x < 2 and < y < 1}. An unforced, undamped autonomous 
version of the double-gyre model was studied by Rom-Kedar, et.al.'^, and an undamped 
system with different forcing was studied by Froyland and Padber^^o] 
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FIG. 2: (Color online) The basin Poincare map for the deterministic part of Eqs. (la)-(lc) at phase 



■0 = 0. The locations of the attracting fixed points corresponding to periodic orbits are denoted 
by the stars. The coloring represents the convergence rate to the attractors. As shown by the 
color bar, positive values converge to the left basin and negative values converge to the right basin. 
The largest magnitude values converge the fastest. Four saddle fixed points are located within 
the boundaries of the domain, and are denoted by large dots. The remaining two are located at 
(0.995,1.005) and (2.010,0.995). 



Prior to defining the control of certain sets of the stochastic system, we first describe 



the important dynamical features of the deterministic part of Eqs. (la)-(lc). There are two 



attracting periodic orbits, which correspond to two fixed points of the Poincare map defined 
by sampling the system at the forcing period. For each fixed point, there corresponds a 
left and a right basin of attraction. Local analysis on the Poincare section about the fixed 
points reveals the attractors to be spiral sinks, which generate the global double-gyre. A 
representative basin map at the phase ip = is shown in Fig. |2} One can see the complicated 
basin boundary structure in which the basins of attraction are intermingled, a signature of 
the existence of a fractal basin boundary. There are also several unstable (saddle) periodic 
orbits corresponding to fixed points that lie along or close to the domain boundary. 

Due to the intermingling of the basin boundaries, one can expect that small perturbations, 
or uncertainties, in initial conditions near the basin boundary will generate large changes 
in dynamical behavior. This can occur deterministically, or when noise is added to the 
system. Therefore, we quantify regions of uncertainty in phase space in both deterministic 
and stochastic settings in Sections IV and V respectively. 
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III. FINITE-TIME LYAPUNOV EXPONENTS 



One method that can be used to understand transport and which quantifies locahzed 
sensitive dependence to initial conditions in a given fluid flow involves the computation of 
finite-time Lyapunov exponents (FTLE). In a deterministic setting, the FTLE also gives 
an explicit measure of phase space uncertainty. Given a dynamical system, one is often 
interested in determining how particles that are initially inflnitesimally close behave as time 
t — )• ±00. It is well-known that a quantitative measure of this asymptotic behavior is 
provided by the classical Lyapunov exponenlpH. In a similar manner, a quantitative measure 
of how much nearby particles separate after a speciflc amount of time has elapsed is provided 
by the FTLE. 

In the early 1990s, PierrehumbertP^l and Pierrehumbert and Yan^^S characterized atmo- 
spheric structures using FTLE flelds. In particular, their work enabled the identiflcation of 
both mixing regions and transport barriers. Later, in a series of papers published in the early 
2000s, Hallei'^^HlS introduced the idea of Lagrangian Coherent Structures (LCS) in order to 
provide a more rigorous, quantitative framework for the identiflcation of fluid structures. 
Hallei'^ proposed that the LCS be deflned as a ridge of the FTLE fleld, and this idea was 
formalized several years later by Shadden, Lekien and Marsdenl^Zl. When computing the 
FTLE fleld of a dynamical system, these LCS, or ridges, are seen to be the structures which 
have a locally maximal FTLE value. 

Although the FTLE/LCS theory can be extended to arbitrary dimensiorpS, in this article 



we consider a 2D velocity fleld : M x / — )• M given by the deterministic part of Eqs. (la) 



(Ic) which is deflned over the time interval / = C M and the following system of 



equations: 

z{t; U, Zo) = v{z{t; U, zq), t), (2a) 
z{ti;ti,Zo) = Zo, (2b) 

where z = {x, yY G M^, Zq G M?, and tel. 

As previously stated, the trajectories of this dynamical system in the inflnite time limits 
can be quantifled with the system's Lyapunov exponents. If one restricts the Lyapunov 
exponent calculation to a flnite time interval, the resulting exponents are the FTLE. In 
practice, the FTLE computation involves consideration of nearby initial conditions and 
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the determination of how the trajectories associated with these initial conditions evolve 
in time. Therefore, the FTLE provides a local measure of sensitivity to initial conditions 
and measures the growth rates of the linearized dynamics about the trajectories. Since the 
details of the derivation of the FTLEp^^ as well as applications that employ the FTLEp^HSS 
have appeared in the literature, we shall only briefly summarize the procedure. 



The solution of the dynamical system given by Eqs. (2a)-(2b) from the initial time ti to 



the final time ti + T can be viewed as the flow map 0*^^^ which is defined as follows: 



M+T . ^ 0*;+^(2o) = z{ti + T; zo). (3) 



We consider an initial point located at z at tj = along with a perturbed point located at 
z + 6z{0) at ti = 0. Using a Taylor series expansion, one finds that 

^^(y) = ^^^^^^(0) + ^(II^^(O)IP)- (4) 
Dropping the higher order terms, the magnitude of the linearized perturbations is given as 

\\6z{T)\\ = y/{6z{0),A), (5) 
where A is the right Cauchy-Green deformation tensor and is given as follows: 

^^ f d4''-^{z{t))Y f d<f)l'-^^{z{t))\ 

with * denoting the adjoint. Then the FTLE can be defined as 

a{z,ti,T) = j^^ln v/A^ax(A), (7) 

where Amax(A) is the maximum eigenvalue of A. 

For a given 2 G at an initial time ti, Eq. ([7| gives the maximum finite-time Lyapunov 
exponent for some finite integration time T (forward or backward), and provides a measure 
of the sensitivity of a trajectory to small perturbations. 

The FTLE field given by a{z,ti,T) can be shown to exhibit ridges of local maxima in 
phase space. The ridges of the field indicate the location of attracting (backward time FTLE 
field) and repelling (forward time FTLE field) structures. In 2D space, the ridge is a curve 
which locally maximizes the FTLE field so that transverse to the ridge one finds the FTLE 
to be a local maximum. 
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FIG. 3: (Color online) Forward FTLE flow field computed using the deterministic form of Eqs. ( la ) 



(Ic) and shown at phase il) = 0. The integration time is T = 20 with an integration step size of 



t = 0.1 and a grid resolution of 0.005 in both x and y. 



Figure |3] shows a snapshot taken at phase = of the forward time FTLE field calculated 



using the deterministic form of Eqs. (la)-(lc) for a finite integration time T = 20. One 



can see in Fig. [3] that there are ridges (in red) of locally maximal FTLE values. These 
ridges or LCS effectively separate the phase space into distinct dynamical regions. For the 
deterministic system, a particle placed in one of these distinct regions will remain in that 
region as the system evolves in time. Notice also that the red ridges appear to be a subset 
of the basin boundary set separating the two basins of attraction shown in Fig. |2} 

In the stochastic system, the noise acts continuously on a particle placed in one of these 
distinct basins, and therefore, it is possible that the particle will cross the LCS and move 
into the other basin. Even though we find the LCS using the deterministic system, the 
location of these structures are a valuable tool in understanding how a particle escapes from 
the basin in which it is initially placed. Coupled with the ideas of uncertain sets and almost 
invariant sets described below, we will use the LCS information provided by the FTLE field 
to determine appropriate regions of control. In this way, it is possible to increase the loitering 
times of a particle in a basin while minimizing the amount of control actuation. 



IV. UNCERTAIN SETS 



In this section, we describe a method that measures the fraction of uncertainty for regions 
in phase space. Nonlinear systems that possess multiple attractors typically can be extremely 
sensitive with respect to the choice of initial conditions. That is, very small changes, or 
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FIG. 4: (Color online) The uncertainty measure for the deterministic part of Eqs. (la)-(lc) for 
the phase ip = 0. A 200 x 100 regular grid of initial conditions is used. The color represents 
the fraction of trajectories that diverge from the trajectory of the center of each ball tested. The 
computation was performed using balls of radius e = 0.02 with twenty points chosen at random 
within each ball. The trajectories were computed until time r = 400 using a threshold parameter 
of 0.8. 



uncertainty, in the initial conditions can lead to different attractors^^"^. Sensitivity here 
is measured in the asymptotic time limit. When applied to high-dimensional attr actors, 
long-time sensitivity is measured with respect to parameter sensitivity'^Sl. In contrast to the 
asymptotic definition, we measure uncertainty in phase space with respect to perturbations 
in initial data by computing exponential changes in distances over short time intervals. 

The sensitivity of the deterministic system to the initial conditions is quantified by defin- 
ing the pointwise uncertainty measure. For a point in phase space (xq, yo) chosen at random, 
a ball of radius e is constructed about that point, B^{xq, yo)- For randomly chosen points 
in the ball B^{xo,yo), we record the number of trajectories, N^{T,xo,yo), at time r, that 
diverge beyond a certain distance from the trajectory of (xo,?/o)- By diverge, we mean that 
we test if the distance between the trajectories at a given time is greater than a threshold 
parameter. The test indicates the fraction of initial conditions in the ball that diverge be- 
yond a certain distance in a given finite time. If there are multiple attractors present, then 
as r — > oo, the uncertain points correspond to those points that go to another basin. For 
finite r, as the radius of the ball is decreased, one obtains more accurate approximations 
of the basin boundary of the system. This is due to the fact that the dynamics are most 
sensitive with respect to the choice of initial conditions when the points straddle the basin 
boundary separating the basins of attraction. 
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Figure |4] shows a sample calculation of the uncertainty measure for the fixed phase if) = 0. 
In this calculation, a 200 x 100 regular grid is used where each grid point coincides with 
the center of a ball of radius 0.02. By comparing the FTLE field shown in Fig. |3]with the 
uncertainty measure shown in Fig. |4| one can see that the location of the FTLE ridge is 
contained within the uncertainty region. In addition, it should be noted that the uncertainty 
measure captures the complicated fractal structure of the basin boundary shown in Fig. |2| 
as radius e approaches zero. For the stochastic double-gyre, it is expected that the regions 
defined by the uncertain sets over all phases, ip G [0,20) will be the most active regions in 
which to actuate if one wishes to increase the residence time of a particle in one particular 
basin. We now refine this notion in the discrete case by defining the actual sets which will 
be almost invariant in the stochastic case. One may think of the almost invariant sets as 
the complement of the most uncertain regions in phase space. 



V. ALMOST INVARIANT SETS 



As mentioned above, approximating the almost invariant sets is another method that can 
quantify where particles loiter in phase space for a stochastic system. Computation of the 
almost invariant sets will be achieved using the stochastic Frobenius- Perron (SFP) operator. 
While the transition probabilities can be approximated for continuous tim^^, we take ad- 
vantage of the fact that the double-gyre system is a periodically forced system. Therefore, 
the dynamics are sampled at a particular phase of the period. For small noise, contributions 
to the probability along trajectories come primarily from the difference between the tangent 
vector and the vector field. Averaging over a period of the drive allows one to sample the 
noise periodically. In doing this, we assume that the noise is not in the tail of the probabil- 
ity distribution so that there are no large, rare events. We therefore construct a Poincare 
map with a specified noise distribution. This method of discrete sampling quantitatively 
identifies the complement of the uncertainty region as two almost invariant sets, which are 
associated with the stable attractors of the associated deterministic system. Although we 
present the machinery for a fixed phase, it is possible to extend the analysis for all phases 
of ip, allowing one to approximate the density over the full forcing period. 

To be explicit, consider a stochastically perturbed map : M" — )• M", x(t) i— i- F{x{t))+r], 
where ?7 is a random variable having the distribution z/(a;). Let F{x{t)) represent the map 
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FIG. 5: (Color online) The results of approximating the almost invariant sets for the stochastic 



system given by Eqs. (la)-(lc) at the phase i/j = 0. A 100 x 100 grid of points is used with 
a standard deviation of o" = 0.05. Figure ^a.) is the mapping of the second eigenvector of the 
reversible matrix back to the phase space. The almost invariant sets are represented by the red 
and the blue regions at the ends of the eigenvector range. Figure [sj^b) is the contour plot of the 
surface, which shows the transition region that forms a barrier between the almost invariant sets. 

on the Poincare section which samples the flow at time intervals of length r, equivalent to 
the period of the steady states. The SFP operator for a normal distribution with mean zero 
and standard deviation a is defined a^ 



Pf[p{x)] = -== / e ^^p{y)dy. (8) 

One approximates the density p(x) by the finite sum of basis functions, p(x) ~ J2iLi '^i4>i{x)i 
where = Xb^x), and xb is an indicator defined on boxes {Bi}f^^ covering the region 

M. The SFP operator is approximated by the N x N matrix, 

Aij = A{B,,B,) = (P^[0,],0i) = / PFM^)]Hx)dx (9) 

J M 



for 1 < i,j < N. Therefore, a transition matrix entry Aij value represents how mass, or 
measure, flows from cell Bi to cell Bj. Details of the method can be found in Ref. [391 

From the transition matrix, a reversible Markov chain is constructed, satisfying vTjAjj = 
T^jAji for all Such a condition implies vr is an invariant (stationary) distribution of 
Since in general A is not reversible as it is defined in Eq. (|9]), a reversible Markov chain is 
constructed from the transition matrix. Let R = {A + A)/2, where Aij = ^Aji and vr is an 
invariant probability density of A. The matrix R now possesses detailed balance. 
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Since R is reversible, it has the properties that its eigenvalues are real and the eigenvectors 
are orthogonal. One can then compute a collection of sets which are almost invariant by 
examining the eigenvectors of R. A gauge of how much measure is ejected each iterate from 
the defined almost invariant sets is given by the eigenvalues. The first few eigenvalues of 
R will be clustered near unity and their associated eigenstates will be a mixture of almost 
invariant sets formed from the above basis elements. Since the first eigenvector of i? is a 
vector of ones, the second eigenvector is used to identify two almost invariant sets emanating 
from basins of attraction. The fact that the eigenvalues cluster near unity means that very 
little mass is ejected from the basins each iterate, which is why the basins are defined as 
almost invariant states. 

The results of the almost invariant sets computation are shown in Fig. |5]for the Poincare 
map of the double- gyre at ip = 0. The second eigenvector can be visualized by a roughly 
piecewise constant function, with a transition region connecting the two pieces. The eigen- 
value associated with this eigenvector is 0.9996. The higher level represents the left basin 
(red) and the lower level represents the right basin (blue). For this computation, a 100 x 100 
grid of points is used with a standard deviation of a = 0.05. As expected, the complement 
set of the left and right almost invariant basins agrees approximately with the location of 
the uncertain region from Fig. |4} 

VI. CORRAL CONTROL IN THE PRESENCE OF FLUCTUATIONS 

Since we have now computed the location of the Lagrangian Coherent Structures and 
have approximated the most uncertain regions and almost invariant regions of phase space, 
we can consider loitering, or residence times, within a defined region. Moreover, we may 
use actuation to control trajectories in the system by increasing the loitering times within 
an almost invariant set. To quantify the residence time of a particle in a basin, we consider 
the double-gyre system with continuously added Gaussian noise. As the standard deviation, 
cr, of the noise is increased, trajectories wander in larger neighborhoods about the stable 
steady states and frequently switch from one basin to another, as shown in Fig. |6] These 
basins are approximated by the almost invariant sets at a fixed phase, shown in Fig. [5j 

The control strategy employs a monitoring ball (that we refer to as a control region) that 
covers the almost invariant set with the center defined as the left non-boundary steady state 
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FIG. 6: (Color online) A graph of the x— coordinate for two simulations of the stochastic system 



given by Eqs. (la)-(lc). For continuous noise with standard deviation a = 0.005, one does not 



observe switching events (red). However, for the larger standard deviation of a = 0.05, one can 
see frequent switching between the two basins (black). 

when time t = 0. When a trajectory passes beyond a threshold distance from the center 
{xsiVs) of the control region, a radial force is turned on. This force moves the trajectory 
back towards the center of the region. This control strategy can be modeled as 



X = —7:Asm{nf{x, t)) cos(7r?/) — ax + rji{t) — c{x — Xs)9(r(x, y) — ro), 

df 

y = TTAcos{TTf{x,t)) sm{TTy) — - ay + ri2{t) - c{y - ys)Q{r{x,y) - ro), 

/(x, t) = e sin(u;t + ip)x^ + (1 — 2e sin(u;t + ■?/'))x, 
r{x,y) = a/(x - + (y _ yJ2^ 



(10a) 

(10b) 

(10c) 
(lOd) 



where c is the magnitude of the control force, tq is the radial threshold, and B is a Heaviside 
function. By adjusting tq and c, it is possible to control the number of actuations and the 
effectiveness of the control scheme. We have designed this specific control strategy to take 
advantage of the underlying flow structure in order to corral the particles into the almost 
invariant set. We therefore refer to our scheme as corral control. 

The control scheme is tested by observing the changes in the almost invariant sets. Fig- 
ure [7] shows the results of approximating the almost invariant sets for the stochastic system 



with control applied to the left basin given by Eqs. (lOa)-(lOd) for c = 0.25. As expected 



there is an increase in the size of the left (red) almost invariant set and a decrease in the 
size of the right (blue) almost invariant set when compared to Fig. [s] 

We have analyzed how individual trajectories escape under this control algorithm. By 
measuring the most frequent path of an escape event from the left basin to the right basin, 
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FIG. 7: (Color online) The results of approximating the almost invariant sets for the stochastic 



system with control applied to the left basin given by Eqs. ( 10a )-( lOd) at the phase ^ = 0. A 100 x 
100 grid of points with a standard deviation of cr = 0.05 is used. Similar to Figure [5]^a), this graph 
is the mapping of the second eigenvector of the reversible matrix back to the phase space. The 
almost invariant sets are represented by the red and the blue regions at the ends of the eigenvector 
range. One can see that the Fig. [Tjleft basin (red) is larger than the Fig. [sj^a) left basin, while the 
Fig. [T] right basin (blue) is smaller than the Fig. [5]^a) right basin. 

one can observe that the trajectories follow the FTLE ridge towards a periodic saddle orbit on 
the X-axis. Topologically, the FTLE ridge loops around the almost invariant sets, corralling 
points in the outlying region to the other side. The goal is to reduce the number of escape 
events by a control algorithm that takes advantage of the topology of the system. 

Figures [8]^a)-(c) illustrate the topology (as does a movie which can be found as auxiliary 
online material). In these figures, 20, 000 particles were placed in the left basin near the fixed 
point and their trajectories were tracked. If a particle crossed the x = 1.2 threshold, the 
particle was considered to have escaped from the left basin into the right basin. In Figs. |8](a)- 
(c), the color contours show the probability density of the trajectory path prehistory for all 
the particles (there are 4361 of them) that escaped from the left basin between phases 
ip G [0.65,0.7). The prehistory shows that there are many paths to escape. However, the 
high probability density in the region of the saddle (near the point (1,0)) shows that most 
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of the particles have been tunneled into the saddle region just prior to escape from the left 
basin. Also shown in the figures are the position of the particles (black dots) before escape, 
the FTLE contour, and a circular region of radius r = 0.4 that denotes the control region. 
These last three features are shown at phase ip = 0.0 in Fig. |8|a), at phase ijj = 0.25 in 
Fig. ^h), and at = 0.5 in Fig. ^c). The control force used in this simulation is c = 0.15. 

In Fig. [sj^a), one can see that although many particles are outside the control region, the 
majority of the particles are inside the control region. For those particles that are inside the 
control region, no control is being applied. In addition, the majority of the particles lie on 
one side of the FTLE ridge, which is bounding the almost invariant set (from comparison 
with Fig. [sj. As time evolves, one can see in Fig. |8](b) that the FTLE ridge sweeps up the 
left side of the figure, and one can see from the prehistory that about half of the particles 
sweep up the left side, staying on the exterior side of the FTLE ridge outside the control 
region so that the control is actuated. The other half of the particles remain on the interior 
side of the FTLE ridge inside the control region so that no control is actuated. 

Eventually, as can be seen in Fig. [sj^c), the FTLE ridge enters the region of uncertainty 
(see Fig. [s]), and all the particles cross over the FTLE ridge. All the particles now have left 
the control region. However, even though the control is being actuated, it cannot overcome 
the underlying dynamics found in this system and all the particles will proceed towards the 
saddle and finally escape. To prevent escape and to increase loitering time, one can use the 
knowledge of the location of the FTLE ridge, uncertainty regions and almost invariant sets 
to design a control region. For example, by choosing a control region with a smaller radius, 
one can prevent the particles from approaching the FTLE ridges and uncertain areas. With 
such a small radius, however, the number of control actuations will be very large. Control 
actuation may be optimized by using a control region with the largest possible radius so 
that the region does not intersect the FTLE structures. 

A. Scaling of the mean escape time with respect to noise and control 

We vary the control parameter, c, to study its effect on the mean escape time. One expects 
that an increase in c will lead to a corresponding increase in the mean time to escape. Note 
that since the control force is finite, there is always a finite probability of escape even if 
control is actuated. The control algorithm decreases the probability of trajectories visiting 
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FIG. 8: (Color online) Probability density of the trajectory path prehistory for particles that 
escaped from the left basin based on a threshold of a; = 1.2 between phases ip G [0.65,0.7). These 
paths form a subset of 20, 000 random particle simulations. Also shown is the position of the 
particles (black dots) before escape, the FTLE ridge (black curve), and the control region of radius 
r = 0.4 (red) at (a) phase ijj = 0.0, (b) phase ip = 0.25, and (c) phase ijj = 0.5. The control force 
is c = 0.15. 

the regions that have a high probability of transition, and the mean time to escape reflects 
the decrease in PDF in the transition regions. 

Figure |9](a) shows the natural log of the average escape time as a function of the inverse 
of the noise intensity. One can see that as the control force is increased, the mean time to 
escape increases with the slope depending linearly on the control parameter. Since small 
changes in c translate into exponential changes in residence times, one recovers a linear 
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FIG. 9: (Color online) (a) Natural log of the average escape time as a function of the inverse of 
the noise intensity. The control radius is r = 0.4 and a threshold of x = 1.2 is used to identify 
escape from the left basin. The results are averaged over 1000 simulations and the slopes of the 
linear fit to each data set is noted on the right, (b) The slopes of the average escape time linear 
fits from Fig. [9]^a) as a function of the control parameter c. 

relation of the slope of the exponent as a function of c which can be seen in Fig. IWb). 



B. Optimizing control actuations 

Another consideration is that one would like to minimize the number of actuations for 
the particle, thereby preserving the energy needed to employ the control scheme. As a 
consequence, this could extend the use of a battery in an autonomous underwater glider. 



Figure [TO] shows the natural log of the average number of actuations per unit time as a 
function of the control radius. For r < 0.36, the FTLE ridge does not intersect the control 
region and the average number of actuations follows the scaling of the average escape time. 
Beyond that, the number of actuations increases as the control scheme works to move the 
trajectory back inside the control region. Therefore, the number of actuations per unit time 
follows the average escape rate for r < 0.36 and the minimum number of actuations can be 
found for the largest control region inside the FTLE ridge. 



C. One-dimensional analysis of control 

To understand the effectiveness of the control scheme and its scaling properties, we ap- 
proximate the natural probability density function (PDF) and the almost invariant set in 
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FIG. 10: Natural log of the average number of actuations per unit time as a function of the control 
radius for parameters a = 0.075 and c = 0.25. The mean number of actuations (o) was calculated 
over 100 simulations of time series with t < 10^, or until an escape event occurred. The error bars 
show one standard deviation of the data from the mean. Overlaid is a linear fit for the data points 
r < 0.36, since for this range of r the control region does not intersect the FTLE ridge. 

the following ID example. Consider a stochastic system, where the deterministic part has 
a stable focus at the origin surrounded by an unstable limit cycle. In the associated deter- 
ministic system, the points inside the limit cycle would spiral toward the stable focus, while 
the points outside would diverge. The stochastic perturbations allow trajectories to escape 
across the limit cycle (and then diverge), creating an almost invariant set. This captures the 
local behavior in one basin of the double-gyre. Now consider a control region with radius 
e > 0, where e determines the distance from the attractor for which the control is actuated. 
Then the system has the form 
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(11) 



where A,a; > 0. As before, c is the magnitude of the control force and G is a Heaviside 

T 

function. Each component of r](t) = [rii(t) , ri2{t)] , is a random variable with intensity D. 

To find the PDF, we switch to polar coordinates using the time-dependent change of 
variables given by x = rcosO and y = rsin6'. This results in the following transformed 
stochastic system: 

r = Ar(r^ ~ 1) ~ crQ{r — e) + Vii't) cosut — r]2(t) smut. (12) 
Here, we find that 6 = —ut and the noise vanishes entirely from the 9 equation. 
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FIG. 11: (Color online) PDF as a function of the control radius r for varying control force c found 



using Eq. (13). Continuous noise is used with a = 1, X = 0.1, and e = 1.5. The normalization 



constant, k, was set in each case so that the area under the curve is one. The inset shows a close-up 
of the PDF near the attractor for the case of no control. 

Assuming the transformed noise term still has intensity D, and k is a normalization 
constant, the probability density p(r), is given by 



p(r) = /cexp ( ^{r\P - 2)) - ^(r^ - e^)Q{r - e) 



(13) 



The PDF is now fully specified by the deterministic flow and noise drift, control strength 
and actuation region, and can be plotted for various parameters. 



In Fig. IT one can see how an increase in the control strength decreases the probability 
of the trajectory moving outside the boundary of the control region. The inset figure shows 
that the local minimum of the PDF occurs at r = 1, which is the location of the unstable 
limit cycle. The PDFs branch at the point of control at r = 1.5. The thick solid curve 
represents the dynamics with no control. It follows that a trajectory escaping outside the 
limit cycle will diverge to infinity, and the PDF increases as one increases r. 

In addition, one can perform the one-dimensional almost invariant set analysis for this 
example. The transport matrix is computed using a grid of 500 intervals on the domain 
and Gaussian noise with standard deviation of cr = 0.1. The point of control is located at 
r = 0.9, to the left of the basin boundary at r = 1. Figure [l2|^a), shows the expansion of 
the left almost invariant set to the right as the control is increased. This information is 
contained in the second eigenvector of the transport matrix. The lower and upper almost 
constant regions of the function represent location of the sets. Figure [l2|(b), shows the 
movement of the transport region to the right as the control is increased. This information 



21 



E 

< 





FIG. 12: (Color online) The results of the almost invariant sets analysis for varying control force 



c in Eq. (12) with A = 0.1, e = 0.9, and Gaussian noise with standard deviation of a = 0.1. 



Figure |l2[a) shows the almost invariant set for varying control force c. The second eigenvector of 
the transport matrix is mapped back to the domain r. Notice that as c increases, the right basin 



extends to the right. Figure 12 b) shows the transition set for varying control force c. The third 
eigenvector of the transport matrix is mapped back to the domain r. Notice that as c increases, 
the transition region moves to the right. 

is contained in the third eigenvector of the transport matrix. The maximum value of the 
function represents the location of the transition region. Notice that it is the complement 
of the almost invariant sets. 

Both the PDF and almost invariant set analysis demonstrate the change that the control 
algorithm has on the natural dynamics of the stochastic system. The control algorithm 
decreases the probability of trajectories visiting regions where deterministic dynamics would 
cause them to diverge. This moves the effective basin boundary farther from the attractor, 
increasing its size. By using sufficient control radius and force, any trajectory that would 
naturally diverge can be redirected towards the attractor. Therefore, the almost invariant 
set can be expanded to the desired size at a cost of the number of control actuations. 

Additionally, it is possible to quantify the relation between the mean escape time and 
the potential defined by the PDF. Using the well-known Kramer's escape ratd^, one can 
predict the rate at which a particle can escape over a potential barrier under Brownian 
motion. In one dimension, the escape time of a particle from a potential defined by U {x) 



IS r 



2tt 



y/U"ix^)\U"{Xf)\ 



exp [{U{xf) — U{xa))/D], where Xa is the location of the attractor and 



X/ is the location of the escape point. From the radial control PDF given by Eq. (13), one 
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can see that the potential function exponent depends hnearly on the control amplitude c. 
All of these one-dimensional analytic results are consistent with the behavior seen for the 
double-gyre system with control that was presented in the preceding subsections. 

VII. CONCLUSIONS 

Using modern dynamical systems theory, we have analyzed the flow of a stochastic double- 
gyre as a model of a wind-driven ocean. The existence of multiple almost invariant sets 
combined with highly uncertain regions in phase space gives rise to stochastic switching 
between the basins. The uncertain regions generate high probability transition sets which 
are the primary cause of switching from one basin to another. We showed that the high 
probability of transitions may be approximated by computing uncertain regions with respect 
to initial conditions and the use of FTLE. 

The system displays maximum sensitivity near the deterministic basin boundary trajec- 
tory. Knowledge of the FTLE ridges in a given flow, in conjunction with knowledge of 
uncertain regions and almost invariant sets, enables one to predict the location of escape 
trajectories from one basin to another. As in previous wor)^, these transition regions may 
be monitored directly to predict future switching events. This knowledge, in turn, allows 
for set-based corral control methods to be used to inhibit the escape event and increase the 
loitering time within the basin. 

Almost invariant sets were computed explicitly using the Frobenius- Perron theory for 
discrete stochastic systems. The sets provided approximate target regions in which to control 
particles which act as sensors in large surveillance regions. That is, given initial conditions in 
one of the almost invariant basins, the particles will stay in the sets for long, but finite time. 
Although computed for a discrete map model, the sets approximate well the continuous 
noise almost invariant sets. They also form sets of the most certain points in that they 
are the complement of the uncertain region, which was computed based on initial condition 
uncertainty. 

Since the almost invariant sets are indicative of long-time, but finite, dynamics, we used 
these sets or subsets to define corral control regions in which we wish to maintain the 
residence times in one basin with minimal control actuation. To this end, we designed a 
simple discrete controller to extend the residence times in the continuous noise double-gyre. 
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We find that even small control amplitudes have an exponential increase on the residence 
times. This is quite similar to the well-known problem of noise- induced escape from a 
potential well, even though the basin boundaries may be fractal in the deterministic case. 

We saw that if the control radius defined a set which did not intersect the ridges defined by 
the FTLE, the number of actuations per unit time followed an exponential law as a function 
of noise intensity. This can be understood completely from the switching rate theory. On 
the other hand, if the control radius intersected the FTLE ridge, the uncertainty regions 
were breached, and the control actuation rate no longer follows a nice scaling law, but rather 
some complicated function of control radius r. 

The model presented here is a simplified version of a double-gyre flow that is a solution 
to a realistic quasi-geostrophic ocean model, and the vehicles controlled are point particles. 
However, the techniques here can be extended to full ocean and glider models in the future. 
The extensions to vehicles with real mass means that inertial effects will need to be included, 
as well extensions to 3D GFD models. However, the machinery presented here is well-suited 
to the description of sets in higher dimensions, and we expect that the monitoring of large 
surveillance regions in the ocean by gliders will be enhanced by implementing our corral 
control method. Additionally, the study and control of coupled systems, with and without 
delay, is of interest. 
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Appendix A: Details of the Ocean Model 

The ocean model results in Fig. [T] were obtained by numerical solution of barotropic 
quasi-geostrophic flow in a single layer basin f2 := [0, 1] x [0, 1] on a /3-plane. The governing 
non-dimensionalized equation for the fluid streamfunction \l/ is: 



+ eJ(^, V^^) + 



(9^ 

dx 



/iV^^ + W, 



(Al) 



dt 
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where J is the Jacobian operator, 

and forcing is provided by a wind stress curl, W , that is prescribed as follows to form a 
double-gyre circulation with a weak periodic "seasonal" variation: 

W = — sin 27iy + 2an sin cut, (A3) 

where the amplitude a = 0.1 and frequency u = 3 were used to produce Fig. [T] 
This system is characterized by the non-dimensional parameter^ 

^=JL '=W^ ^^^^ 

where /3 is the rotation parameter, R is the bottom friction and L and U are respectively the 
characteristic length and velocity scales of the basin. The parameters /i and e correspond to 
the relative length scales of the Stommel and inertial layers, respectively. To produce Fig. [TJ 
we used fi = 0.04 and e = 0.0004. 

The above model has been numerically integrated using second-order spatial differences 
and second-order Runge-Kutta time stepping for the streamfunction, and a fourth-order 
Runge-Kutta algorithm to compute the Lagrangian trajectories of tracer particles. In Fig. [l| 
a grid of resolution 64^ was used and a dimensionless time step of dt = 0.001. Coordinates 
of tracer particles are independent of the grid while flow velocities at the particle locations 
are found using bilinear interpolation from the grid values. An initially static ocean spins-up 
in a few hundred time steps. If a = the spun-up solution is stationary, while non-zero a 
leads to a superimposed oscillatory behavior. The tracers are held in place until spin-up is 
complete. 
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